
***************************************************************
*  Estimating LLF effects on any sex preference strategies, 
*  stopping rule use and postnatal selection
***************************************************************	

***************************************************
* (1.0) Import and prep the data
***************************************************	
	
	use 	"$f_2PKbirthpanel", clear	
	
	do `"$d_do\LLF_StackData.do"'
	
	cap drop years_since_LLF
	gen years_since_LLF = year-LLFyear if LLFyear != .
	
	rename order_livebirth order 

	gen missingmale = 1 if cmale==.
	bysort id: egen hasmissing = total(missingmale)
	drop if hasmissing !=0

	drop if ultrasound == 1
	
		cap drop period 
		gen period = . 
		replace period = 1 if year>=1964 & years_since_LLF <= 0
		replace period = 2 if years_since_LLF >= 1 & years_since_LLF <= 4
		replace period = 3 if years_since_LLF >= 5 & year<=1979 & years_since_LLF !=.
		
		cap drop period_* 
		tab period, gen(period_)
	
		cap drop period2_male period3_male
		gen period2_male = period_2*cmale
		gen period3_male = period_3*cmale	
		
		cap drop period2_noson period3_noson 
		gen period2_noson = period_2*noson
		gen period3_noson = period_3*noson
	
		gen cmale_noson = cmale*noson
		gen period2_male_noson = period2_male*noson
		gen period3_male_noson = period3_male*noson

	cap drop terminal 
	gen terminal = 0
	replace terminal = 1 if order == tot_kids

	bysort id: egen temp1 = mean(cmale) if order==1
	bysort id: egen temp2 = mean(cmale) if order==2
	bysort id: egen temp3 = mean(cmale) if order==3
	bysort id: egen temp4 = mean(cmale) if order==4
	bysort id: egen temp5 = mean(cmale) if order==5
	bysort id: egen temp6 = mean(cmale) if order==6
	
	bysort id: egen first_male = mean(temp1) 
	bysort id: egen second_male = mean(temp2) 
	bysort id: egen third_male = mean(temp3) 
	bysort id: egen fourth_male = mean(temp4) 
	bysort id: egen fifth_male = mean(temp5) 
	bysort id: egen sixth_male = mean(temp6) 

	replace han = 0 if han == .

	merge m:1 prov year using  "$d_data\IMR.dta", gen(IMRmerge)

	cap drop _merge
	merge 	m:1 prov year using `"$d_data\ProductionData.dta"'
	
	cap drop _merge
	merge m:1 prov using `"$d_data\sentdown_mean.dta"'
	drop _merge 
	
	gen 	culturalrev = 0 
	replace culturalrev = mean_sentdown if year >= 1966 & year <=1969

	global controls1 	"age_marriage i.educ U5MR_5yr_ave"
	global controls2 	"grain gdp pop_total_agric agric_production culturalrev"
	
	forval i = 1/5 { 
		cap drop surv_`i' 
		gen surv_`i' = 1 if cmale != . 
			replace surv_`i' = 0 if child_death_age < `i' & cmale != .
		} 	
	
	drop if year <= 1964 
	drop if year >= 1980
	drop if year == 1968
	drop if year == 1967
			
	****************************
	replace order = 3 if order>3 & order!=. 
	
	gen order_2 = order==2
	gen order_3 = order==3
	gen order2_noson = order_2*noson
	gen order3_noson = order_3*noson
	
	drop if years_since_LLF>7
	char FE_prov[omit] 108 
	drop if years_since_LLF<-8

	keep if rural==1

	set seed 8675309

*****************************************************
* (2.0) Sex selection regressions 
*****************************************************
	
	***************************************************
	* (2.1) First Birth: Effect of the sex of the
	* 		birth on the probability of stopping
	***************************************************	

		*************************************************
		* Fertility Discontinuation 
		* (probability a given birth is the last birth)
		*************************************************

		xi: reg terminal cmale period_2 period_3 period2_male period3_male ///
			$controls1 $controls2 i.FE_t i.FE_prov ///
			if year!=1967 & year!=1968 & order==1, cluster(prov)
			
		* Hypothesis tests: we are interested in knowing if the sex of the birth 
		* affects the probability of stopping - above and beyond the baseline
		* probability of stopping in any period (ie allowing for the probability
		* of stopping to change over time). So sex of the child + sex*period interactions
		
				* Period: Pre LLF
				lincom cmale + 0 
				local est_p1 = r(estimate)
				
				boottest (cmale = 0), bootcluster(prov) weighttype(rademacher) ///
				boottype(wild) reps(999) nograph gridpoints(100)
				mat ci = r(CI)
				local pval =  r(p)

				if `pval' < 0.1 { 
					local pre_noson = `est_p1' 
					}
				else if `pval' >= 0.1 { 
					local pre_noson = 0  
					}	
				
				* Store results for Table 1
				mat T1_noson = `est_p1', . \ `pval' , . \  ci[1,1], ci[1,2] 
				mat T1_son = J(3,6,.)
					
				* Period: Early LLF
				lincom cmale + period2_male 
				local est_e1 = r(estimate)
				
				boottest (cmale + period2_male = 0), bootcluster(prov) weighttype(rademacher) ///
				boottype(wild) reps(999) nograph gridpoints(100)
				mat ci = r(CI)
				local pval =  r(p)
				
				if `pval' < 0.1 { 
					local early_noson = `est_e1' 
					}
				else if `pval' >= 0.1 { 
					local early_noson = 0  
					}	
					
				mat temp = `est_e1', . \ `pval' , . \  ci[1,1], ci[1,2] 
				mat T1_noson = T1_noson, temp 
				
				* Period: Late LLF
				lincom cmale + period3_male 
				local est_l1 = r(estimate)
				
				boottest (cmale + period3_male = 0), bootcluster(prov) weighttype(rademacher) ///
				boottype(wild) reps(999) nograph gridpoints(100)	
				mat ci = r(CI)
				local pval =  r(p)
							
				if `pval' < 0.1 { 
					local late_noson = `est_l1'  
					}
				else if `pval' >= 0.1 { 
					local late_noson = 0  
					}	
					
				mat temp = `est_l1', . \ `pval' , . \  ci[1,1], ci[1,2]
				mat T1_noson = T1_noson, temp 
				
			* Save the point estimates here for population-level incidence calculations later
			* treating statistically insignificant point estimates as zero
			* Col 2 contains estimates for couples with no previously born son, col 2 for couples with at least one son
			* Rows correspond to LLF periods (pre, early, late)
			mat stop_pre = `pre_noson', 0
			mat stop_early = `early_noson', 0
			mat stop_late = `late_noson', 0 
			
			* Count of births and terminal births in each period, needed to compute weights
			* for population-level incidence calculations later 
			count if period == 1 & e(sample) 
				local pre_noson = r(N) 
			count if period == 2 & e(sample) 
				local early_noson = r(N)
			count if period == 3 & e(sample) 
				local late_noson = r(N)
							
				mat n_pre = `pre_noson', 0
				mat n_early = `early_noson', 0
				mat n_late = `late_noson', 0 

			count if period == 1 & terminal == 1 & e(sample) 
				local pre_noson = r(N) 
			count if period == 2 & terminal == 1 & e(sample) 
				local early_noson = r(N)
			count if period == 3 & terminal == 1 & e(sample) 
				local late_noson = r(N)
				
				mat c_pre = `pre_noson', 0
				mat c_early = `early_noson', 0
				mat c_late = `late_noson', 0 
					
		*******************************
		* Postnatal Selection 
		* - overall probability a given birth is male
		*******************************
		
		xi: reg cmale period_2 period_3  ///
			$controls1 $controls2 i.FE_t i.FE_prov ///
			if year!=1967 & year!=1968 & order==1, cluster(prov)
			estimates store preferred_1
			
			*Couples with sons not estimated because these are first parity and no couples have a son
			*For this reason there is not a comparison group available to us that 
			*would permit estimates comparable to higher parity birth estimates we report.
			*These are effect of period on likelihood of a male birth relative to pre-LLF first parity births. 
			
			*Instead, we focus on the effect of not having a son on the probaiblity of a male births
			*relative to births to couples with a son in the same parity, same period. 
			
			*We don't report these estimates, but we have to set the values to zero so that we can do 
			* incidence calculations later on. 
			
			mat postnatal_pre = 0, 0
			mat postnatal_early = 0, 0
			mat postnatal_late = 0, 0

		* Controlling for total number of kids 
			xi: reg cmale period_2 period_3  ///
			$controls1 $controls2 tot_kids i.FE_t i.FE_prov ///
			if year!=1967 & year!=1968 & order==1, cluster(prov)
			estimates store robust_1
			
	************************
	* (2.2) Second Birth
	************************

		*************************************************
		*  Fertility Discontinuation 
		* (probability a given birth is the last birth)
		*************************************************
		
		xi: reg terminal cmale noson period_2 period_3  ///
			cmale_noson period2_noson period3_noson ///
			period2_male period3_male period2_male_noson period3_male_noson ///
			$controls1 $controls2 i.FE_t i.FE_prov if year!=1967 & year!=1968 & order==2, cluster(prov)
	
		* Hypothesis tests: we are interested in knowing if the sex of the birth 
		* affects the probability of stopping
		
				* Period: Pre LLF; Previous Son: No 	 
				lincom cmale + cmale_noson
				local est_p2_ns = r(estimate)
				
				boottest (cmale + cmale_noson = 0), bootcluster(prov) weighttype(rademacher) ///
				boottype(wild) reps(999) nograph gridpoints(100)
				mat ci = r(CI)
				local pval =  r(p)
			
				if `pval' < 0.1 { 
					local pre_noson = `est_p2_ns' 
					}
				else if `pval' >= 0.1 { 
					local pre_noson = 0  
					}	
	
				mat temp1 = `est_p2_ns', . \ `pval', . \  ci[1,1], ci[1,2] 
	
				* Period: Pre LLF; Previous Son: Yes	 
				lincom cmale + 0 
				local est_p2_s = r(estimate)
				
				boottest (cmale = 0), bootcluster(prov) weighttype(rademacher) ///
				boottype(wild) reps(999) nograph gridpoints(100)
				
				mat ci = r(CI)
				local pval =  r(p)

				if `pval' < 0.1 { 
					local pre_son =`est_p2_s' 
					}
				else if `pval' >= 0.1 { 
					local pre_son = 0  
					}
					
				mat temp2 = `est_p2_s', . \ `pval', . \  ci[1,1], ci[1,2] 
		
				* Period: Early LLF; Previous Son: No 	 
				lincom  cmale + cmale_noson + period2_male + period2_male_noson
				local est_e2_ns = r(estimate)

				boottest (cmale + cmale_noson + period2_male + period2_male_noson = 0), ///
				bootcluster(prov) weighttype(rademacher) ///
				boottype(wild) reps(999) nograph gridpoints(100)
				
				mat ci = r(CI)
				local pval =  r(p)
			
				if `pval' < 0.1 { 
					local early_noson = `est_e2_ns' 
					}
				else if `pval' >= 0.1 { 
					local early_noson = 0  
					}
					
				mat temp = `est_e2_ns', . \ `pval', . \  ci[1,1], ci[1,2]
				mat temp1 = temp1, temp

				* Period: Early LLF; Previous Son: Yes 	 
				lincom cmale + period2_male 
				local est_e2_s = r(estimate)
				
				boottest (cmale + period2_male = 0), bootcluster(prov) weighttype(rademacher) ///
				boottype(wild) reps(999) nograph gridpoints(100)
				
				mat ci = r(CI)
				local pval =  r(p)
				
				if `pval' < 0.1 { 
					local early_son = `est_e2_s'  
					}
				else if `pval' >= 0.1 { 
					local early_son = 0  
					}	
					
				mat temp = `est_e2_s', . \ `pval', . \  ci[1,1], ci[1,2]
				mat temp2 = temp2, temp
					
				* Period: Late LLF; Previous Son: No 	 
				lincom cmale + cmale_noson + period3_male + period3_male_noson
				local est_l2_ns = r(estimate)

				boottest (cmale + cmale_noson + period3_male + period3_male_noson = 0), ///
				bootcluster(prov) weighttype(rademacher) ///
				boottype(wild) reps(999) nograph gridpoints(100)
				mat ci = r(CI)
				local pval =  r(p)
			
				if `pval' < 0.1 { 
					local late_noson = `est_l2_ns'  
					}
				else if `pval' >= 0.1 { 
					local late_noson = 0  
					}

				mat temp = `est_l2_ns', . \ `pval', . \  ci[1,1], ci[1,2]
				mat temp1 = temp1, temp
					
				* Period: Late LLF; Previous Son: Yes 	 
				lincom cmale + period3_male 
				local est_l2_s = r(estimate)
				
				boottest (cmale + period3_male = 0), bootcluster(prov) weighttype(rademacher) ///
				boottype(wild) reps(999) nograph gridpoints(100)
				mat ci = r(CI)
				local pval =  r(p)
			
				if `pval' < 0.1 { 
					local late_son = `est_l2_s'  
					}
				else if `pval' >= 0.1 { 
					local late_son = 0  
					}

				mat temp = `est_l2_s', . \ `pval', . \  ci[1,1], ci[1,2]
				mat temp2 = temp2, temp
			
			mat T1_noson = T1_noson \ temp1
			mat T1_son = T1_son \ temp2
			
			* Store the point estimates for each group in each period
			mat temp1 = `pre_noson', `pre_son'
			mat temp2 = `early_noson', `early_son'
			mat temp3 = `late_noson', `late_son' 

			mat stop_pre = stop_pre \ temp1
			mat stop_early = stop_early \ temp2
			mat stop_late = stop_late \ temp3
		
			* Store the count of births and terminal births in each period 
			count if period == 1 & noson == 1 & e(sample) 
				local pre_noson = r(N) 
			count if period == 2 & noson == 1 & e(sample) 
				local early_noson = r(N)
			count if period == 3 & noson == 1 & e(sample) 
				local late_noson = r(N)
				
			count if period == 1 & noson == 0 & e(sample) 
				local pre_son = r(N) 
			count if period == 2 & noson == 0 & e(sample) 
				local early_son = r(N)
			count if period == 3 & noson == 0 & e(sample) 
				local late_son = r(N)
				
				mat temp1 = `pre_noson', `pre_son'
				mat temp2 = `early_noson', `early_son'
				mat temp3 = `late_noson', `late_son' 

				mat n_pre = n_pre \ temp1
				mat n_early = n_early \ temp2
				mat n_late = n_late \ temp3

			count if period == 1 & noson == 1 & terminal == 1 & e(sample) 
				local pre_noson = r(N) 
			count if period == 2 & noson == 1 & terminal == 1 & e(sample) 
				local early_noson = r(N)
			count if period == 3 & noson == 1 & terminal == 1 & e(sample) 
				local late_noson = r(N)
				
			count if period == 1 & noson == 0 & terminal == 1 & e(sample) 
				local pre_son = r(N) 
			count if period == 2 & noson == 0 & terminal == 1 & e(sample) 
				local early_son = r(N)
			count if period == 3 & noson == 0 & terminal == 1 & e(sample) 
				local late_son = r(N)
				
				mat temp1 = `pre_noson', `pre_son'
				mat temp2 = `early_noson', `early_son'
				mat temp3 = `late_noson', `late_son' 

				mat c_pre = c_pre \ temp1
				mat c_early = c_early \ temp2
				mat c_late = c_late \ temp3
						
		*******************************
		* Postnatal Selection
		*******************************

		xi: reg cmale noson period_2 period_3 period2_noson period3_noson   ///
			$controls1 $controls2 i.FE_t i.FE_prov if ///
			year!=1967 & year!=1968 & order == 2, cluster(prov)
			estimates store preferred_2

			* Hypothesis tests: Interested in the effect of the period and 
			* lack of a son on probability of having a son, relative to births 
			* in the same period at same parity
				* Note births with a son in the same period is the reference group
				* Period: Pre LLF; Previous son: No
				lincom noson + 0
				local est_p2_ns = r(estimate)
				
				boottest (noson = 0), bootcluster(prov) weighttype(rademacher) ///
				boottype(wild) reps(999) nograph gridpoints(100)

				mat ci = r(CI)
				local pval =  r(p)
				
				if `pval' < 0.1 { 
					local pre_noson = `est_p2_ns' 
					}
				else if `pval' >= 0.1 { 
					local pre_noson = 0  
					}	

				mat T2 = `est_p2_ns', . \ `pval', . \  ci[1,1], ci[1,2] 
				mat TA9_p1 = `est_p2_ns', . \ `pval', . \ ci[1,1], ci[1,2] 
				
				* Period: Early LLF; Previous son: No
				lincom noson + period2_noson
				local est_e2_ns = r(estimate)
				
				boottest (noson + period2_noson = 0), bootcluster(prov) weighttype(rademacher) ///
				boottype(wild) reps(999) nograph gridpoints(100)

				mat ci = r(CI)
				local pval =  r(p)
			
				if `pval' < 0.1 { 
					local early_noson = `est_e2_ns' 
					}
				else if `pval' >= 0.1 { 
					local early_noson = 0  
					}	

				mat temp = `est_e2_ns', . \ `pval', . \  ci[1,1], ci[1,2] 
				mat T2 = T2, temp 
				mat TA9_p1 = TA9_p1 \ `est_e2_ns', . \ `pval', . \  ci[1,1], ci[1,2] 
				
				* Period: Late LLF; Previous son: No
				lincom noson + period3_noson
				local est_l2_ns = r(estimate)
				
				boottest (noson + period3_noson = 0), bootcluster(prov) weighttype(rademacher) ///
				boottype(wild) reps(999) nograph gridpoints(100)

				mat ci = r(CI)
				local pval =  r(p)
			
				if `pval' < 0.1 { 
					local late_noson = `est_l2_ns'  
					}
				else if `pval' >= 0.1 { 
					local late_noson = 0 
					}	

				mat temp = `est_l2_ns', . \ `pval', . \  ci[1,1], ci[1,2]
				mat T2 = T2, temp 
				mat TA9_p1 = TA9_p1 \ `est_l2_ns', . \ `pval', . \  ci[1,1], ci[1,2] 
					
			mat temp1 = `pre_noson', 0 
			mat temp2 = `early_noson', 0 
			mat temp3 = `late_noson', 0 

			mat postnatal_pre = postnatal_pre \ temp1
			mat postnatal_early = postnatal_early \ temp2
			mat postnatal_late = postnatal_late \ temp3

		*******************************
		* Postnatal Neglect 
		* Survival to ages 1-5
		*******************************
		
		forval age = 1/5 { 
			di "************************************"
			di "*******Second Births Age `age'********"
			di "************************************"
			
			xi: reg cmale noson period_2 period_3 period2_noson period3_noson   ///
				$controls1 $controls2 i.FE_t i.FE_prov if year!=1967 & year!=1968 ///
				& order == 2 & surv_`age'==1, cluster(prov)

				lincom noson + 0
				local est_p = r(estimate)
					boottest (noson = 0), bootcluster(prov) weighttype(rademacher) ///
					boottype(wild) reps(999) nograph gridpoints(100)
					mat ci_p = r(CI)
					local pval_p =  r(p)

				lincom noson + period2_noson
				local est_e = r(estimate)
					boottest (noson + period2_noson = 0), bootcluster(prov) weighttype(rademacher) ///
					boottype(wild) reps(999) nograph gridpoints(100)
					mat ci_e = r(CI)
					local pval_e =  r(p)

				lincom noson + period3_noson
				local est_l = r(estimate)
					boottest (noson + period3_noson = 0), bootcluster(prov) weighttype(rademacher) ///
					boottype(wild) reps(999) nograph gridpoints(100)
					mat ci_l = r(CI)
					local pval_l =  r(p)
					
				mat temp1 = `est_p', . \ `pval_p', . \ ci_p[1,1], ci_p[1,2] 
				mat temp2 = `est_e', . \ `pval_e', . \ ci_e[1,1], ci_e[1,2] 
				mat temp3 = `est_l', . \ `pval_l', . \ ci_l[1,1], ci_l[1,2]
				mat temp =  temp1 \ temp2 \ temp3 
				mat TA9_p1 = TA9_p1, temp
			}

			
		* Controlling for total number of kids 
		xi: reg cmale noson period_2 period_3 period2_noson period3_noson   ///
			$controls1 $controls2 tot_kids i.FE_t i.FE_prov if year!=1967 ///
			& year!=1968 & order == 2 , cluster(prov)
			estimates store robust_2
		
	************************
	* (2.3) Third + Births
	************************

		*************************************************
		*  Fertility Discontinuation 
		* (probability a given birth is the last birth)
		*************************************************

		xi: reg terminal cmale noson period_2 period_3  ///
			cmale_noson period2_noson period3_noson period2_male period3_male ///
			period2_male_noson period3_male_noson $controls1 $controls2 ///
			i.FE_t i.FE_prov if year!=1967 & year!=1968 & order>=3, cluster(prov)
		
		* Hypothesis tests: we are interested in knowing if the sex of the birth 
		* affects the probability of stopping
				* Period: Pre LLF; Previous Son: No 	 
					lincom cmale + cmale_noson
					local est_p3_ns = r(estimate)

					boottest (cmale + cmale_noson = 0), bootcluster(prov) weighttype(rademacher) ///
					boottype(wild) reps(999) nograph gridpoints(100)
					
					mat ci = r(CI)
					local pval =  r(p)				
				
					if `pval' < 0.1 { 
						local pre_noson = `est_p3_ns'  
						}
					else if `pval' >= 0.1 { 
						local pre_noson = 0  
						}	
						
				mat temp1 = `est_p3_ns', . \ `pval', . \  ci[1,1], ci[1,2] 
					
				* Period: Pre LLF; Previous Son: Yes 	 
					lincom cmale + 0 
					local est_p3_s = r(estimate)
					
					boottest (cmale = 0), bootcluster(prov) weighttype(rademacher) ///
					boottype(wild) reps(999) nograph gridpoints(100)
			
					mat ci = r(CI)
					local pval =  r(p)				
				
					if `pval' < 0.1 { 
						local pre_son = `est_p3_s'  
						}
					else if `pval' >= 0.1 { 
						local pre_son = 0  
						}	

				mat temp2 = `est_p3_s', . \ `pval', . \  ci[1,1], ci[1,2] 
						
				* Period: Early LLF; Previous Son: No 	 
					lincom  cmale + cmale_noson + period2_male + period2_male_noson
					local est_e3_ns = r(estimate)
					
					boottest (cmale + cmale_noson + period2_male + period2_male_noson = 0), ///
					bootcluster(prov) weighttype(rademacher) ///
					boottype(wild) reps(999) nograph gridpoints(100)

					mat ci = r(CI)
					local pval =  r(p)				
				
					if `pval' < 0.1 { 
						local early_noson = `est_e3_ns'  
						}
					else if `pval' >= 0.1 { 
						local early_noson = 0  
						}
	
				mat temp = `est_e3_ns', . \ `pval', . \  ci[1,1], ci[1,2] 
				mat temp1 = temp1, temp
				
				* Period: Early LLF; Previous Son: Yes 	 
					lincom cmale + period2_male 
					local est_e3_s = r(estimate)
					
					boottest (cmale + period2_male = 0), bootcluster(prov) weighttype(rademacher) ///
					boottype(wild) reps(999) nograph gridpoints(100)

					mat ci = r(CI)
					local pval =  r(p)				
				
					if `pval' < 0.1 { 
						local early_son = `est_e3_s'  
						}
					else if `pval' >= 0.1 { 
						local early_son = 0  
						}	
						
				mat temp = `est_e3_s', . \ `pval', . \  ci[1,1], ci[1,2] 
				mat temp2 = temp2, temp
						
				* Period: Late LLF; Previous Son: No 	 
					lincom cmale + cmale_noson + period3_male + period3_male_noson
					local est_l3_ns = r(estimate)
					
					boottest (cmale + cmale_noson + period3_male + period3_male_noson = 0), ///
					bootcluster(prov) weighttype(rademacher) boottype(wild) reps(999) ///
					nograph gridpoints(100)

					mat ci = r(CI)
					local pval =  r(p)				
				
					if `pval' < 0.1 { 
						local late_noson = `est_l3_ns'  
						}
					else if `pval' >= 0.1 { 
						local late_noson = 0  
						}

				mat temp = `est_l3_ns', . \ `pval', . \  ci[1,1], ci[1,2]
				mat temp1 = temp1, temp
						
				* Period: Late LLF; Previous Son: Yes 	 
					lincom cmale + period3_male 
					local est_l3_s = r(estimate)
					
					boottest (cmale + period3_male = 0), bootcluster(prov) weighttype(rademacher) ///
					boottype(wild) reps(999) nograph gridpoints(100)
		
					mat ci = r(CI)
					local pval =  r(p)				
				
					if `pval' < 0.1 { 
						local late_son = `est_l3_s' 
						}
					else if `pval' >= 0.1 { 
						local late_son = 0  
						}
						
				mat temp = `est_l3_s', . \ `pval', . \  ci[1,1], ci[1,2]
				mat temp2 = temp2, temp
				
				mat T1_noson = T1_noson \ temp1
				mat T1_son = T1_son \ temp2
						
				* Store the point estimates by previous son, for each period
				mat temp1 = `pre_noson', `pre_son'
				mat temp2 = `early_noson', `early_son'
				mat temp3 = `late_noson', `late_son' 

				mat stop_pre = stop_pre \ temp1
				mat stop_early = stop_early \ temp2
				mat stop_late = stop_late \ temp3
					
				* Count of births and terminal births in each period 
				count if period == 1 & noson == 1 & e(sample) 
					local pre_noson = r(N) 
				count if period == 2 & noson == 1 & e(sample) 
					local early_noson = r(N)
				count if period == 3 & noson == 1 & e(sample) 
					local late_noson = r(N)
					
				count if period == 1 & noson == 0 & e(sample) 
					local pre_son = r(N) 
				count if period == 2 & noson == 0 & e(sample) 
					local early_son = r(N)
				count if period == 3 & noson == 0 & e(sample) 
					local late_son = r(N)
					
					mat temp1 = `pre_noson', `pre_son'
					mat temp2 = `early_noson', `early_son'
					mat temp3 = `late_noson', `late_son' 

					mat n_pre = n_pre \ temp1
					mat n_early = n_early \ temp2
					mat n_late = n_late \ temp3

				count if period == 1 & noson == 1 & terminal == 1 & e(sample) 
					local pre_noson = r(N) 
				count if period == 2 & noson == 1 & terminal == 1 & e(sample) 
					local early_noson = r(N)
				count if period == 3 & noson == 1 & terminal == 1 & e(sample) 
					local late_noson = r(N)
					
				count if period == 1 & noson == 0 & terminal == 1 & e(sample) 
					local pre_son = r(N) 
				count if period == 2 & noson == 0 & terminal == 1 & e(sample) 
					local early_son = r(N)
				count if period == 3 & noson == 0 & terminal == 1 & e(sample) 
					local late_son = r(N)
					
					mat temp1 = `pre_noson', `pre_son'
					mat temp2 = `early_noson', `early_son'
					mat temp3 = `late_noson', `late_son' 

					mat c_pre = c_pre \ temp1
					mat c_early = c_early \ temp2
					mat c_late = c_late \ temp3
					
			*******************************
			* Postnatal Selection
			*******************************
					
			xi: reg cmale noson period_2 period_3 period2_noson period3_noson   ///
				$controls1 $controls2 i.FE_t i.FE_prov if year!=1967 & year!=1968 ///
				& order >= 3 & years_since_LLF<8, cluster(prov)
				estimates store preferred_3
				
			* Hypothesis tests: Interested in the effect of the period and 
			* lack of a son on probability of having a son, relative to births 
			* in the same period at same parity
				* Note births with a son in the same period is the reference group

				* Period: Pre LLF; Previous son: No
				lincom noson + 0
				local est_p3_ns = r(estimate)
				
				boottest (noson = 0), bootcluster(prov) weighttype(rademacher) ///
				boottype(wild) reps(999) nograph gridpoints(100)

				mat ci = r(CI)
				local pval =  r(p)
			
				if `pval' < 0.1 { 
					local pre_noson = `est_p3_ns'  
					}
				else if `pval' >= 0.1 { 
					local pre_noson = 0  
					}	
					
				mat temp1 = `est_p3_ns', . \ `pval', . \  ci[1,1], ci[1,2] 
				mat TA9_p2 = `est_p3_ns', . \ `pval', . \ ci[1,1], ci[1,2] 
				
				* Period: Early LLF; Previous son: No
				lincom noson + period2_noson
				local est_e3_ns = r(estimate)

				boottest (noson + period2_noson = 0), bootcluster(prov) weighttype(rademacher) ///
				boottype(wild) reps(999) nograph gridpoints(100)
				
				mat ci = r(CI)
				local pval =  r(p)
			
				if `pval' < 0.1 { 
					local early_noson = `est_e3_ns'  
					}
				else if `pval' >= 0.1 { 
					local early_noson = 0  
					}	

				mat temp = `est_e3_ns', . \ `pval', . \  ci[1,1], ci[1,2] 
				mat temp1 = temp1, temp
				mat TA9_p2 = TA9_p2 \ `est_e3_ns', . \ `pval', . \ ci[1,1], ci[1,2] 
				
				* Period: Late LLF; Previous son: No
				lincom noson + period3_noson
				local est_l3_ns = r(estimate)
	
				boottest (noson + period3_noson = 0), bootcluster(prov) weighttype(rademacher) ///
				boottype(wild) reps(999) nograph gridpoints(100)
	
				mat ci = r(CI)
				local pval =  r(p)
			
				if `pval' < 0.1 { 
					local late_noson = `est_l3_ns'   
					}
				else if `pval' >= 0.1 { 
					local late_noson = 0 
					}	

				mat temp = `est_l3_ns', . \ `pval', . \  ci[1,1], ci[1,2] 
				mat temp1 = temp1, temp
				mat TA9_p2 = TA9_p2 \ `est_l3_ns', . \ `pval', . \ ci[1,1], ci[1,2] 
			
			mat T2 = T2 \ temp1 
			
			mat temp1 = `pre_noson', 0 
			mat temp2 = `early_noson', 0
			mat temp3 = `late_noson', 0 

			mat postnatal_pre = postnatal_pre \ temp1
			mat postnatal_early = postnatal_early \ temp2
			mat postnatal_late = postnatal_late \ temp3		

		*******************************
		* Postnatal Neglect 
		* Survival to ages 1-5
		*******************************
		forval age = 1/5 { 
			di "************************************"
			di "*******Third Births Age `age'********"
			di "************************************"
			
			xi: reg cmale noson period_2 period_3 period2_noson period3_noson   ///
				$controls1 $controls2 i.FE_t i.FE_prov if year!=1967 & year!=1968 ///
				& order == 3 & surv_`age'==1, cluster(prov)

				lincom noson + 0
				local est_p = r(estimate)
					boottest (noson = 0), bootcluster(prov) weighttype(rademacher) ///
					boottype(wild) reps(999) nograph gridpoints(100)
					mat ci_p = r(CI)
					local pval_p =  r(p)

				lincom noson + period2_noson
				local est_e = r(estimate)
					boottest (noson + period2_noson = 0), bootcluster(prov) weighttype(rademacher) ///
					boottype(wild) reps(999) nograph gridpoints(100)
					mat ci_e = r(CI)
					local pval_e =  r(p)

				lincom noson + period3_noson
				local est_l = r(estimate)
					boottest (noson + period3_noson = 0), bootcluster(prov) weighttype(rademacher) ///
					boottype(wild) reps(999) nograph gridpoints(100)
					mat ci_l = r(CI)
					local pval_l =  r(p)
					
				mat temp1 = `est_p', . \ `pval_p', . \ ci_p[1,1], ci_p[1,2] 
				mat temp2 = `est_e', . \ `pval_e', . \ ci_e[1,1], ci_e[1,2] 
				mat temp3 = `est_l', . \ `pval_l', . \ ci_l[1,1], ci_l[1,2]
				mat temp =  temp1 \ temp2 \ temp3 
				mat TA9_p2 = TA9_p2, temp
		}
				
		* Controlling for total number of kids 
		xi: reg cmale noson period_2 period_3 period2_noson period3_noson   ///
			$controls1 $controls2 tot_kids i.FE_t i.FE_prov if year!=1967 ///
			& year!=1968 & order >= 3 & years_since_LLF<8, cluster(prov)
			estimates store robust_3
		
			
********************************************
* (3.0) Format and output tables 
********************************************
		
	********************************************
	* (3.1) Table 1
	********************************************
		mat T1 = T1_noson, T1_son
		putexcel set  `"$d_tab\Table_1.xlsx"', replace 	
		local excelcols "B C D F G H"
		
		local xc = 1
		forval c = 1(2)12 { 
			local xr = 3
			forval r = 1(3)9 {
				local c2 = `c'+1
				local r2 = `r'+1
				local r3 = `r'+2
				if T1[`r2',`c']<.1 {
					local lincom = round(T1[`r',`c'],.001)
					local lincom = string(`lincom', "%05.3f")
					scalar pt_est = "`lincom'**"
					}
				if T1[`r2',`c']<.05 {
					local lincom = round(T1[`r',`c'],.001)
					local lincom = string(`lincom', "%05.3f")
					scalar pt_est = "`lincom'**"
					}
				if T1[`r2',`c']<.001 {
					local lincom = round(T1[`r',`c'],.001)
					local lincom = string(`lincom', "%05.3f")
					scalar pt_est = "`lincom'***"
					}
				else if T1[`r2',`c']>=.1 {
					local lincom = round(T1[`r',`c'],.001)
					local lincom = string(`lincom', "%05.3f")
					scalar pt_est = "`lincom'"
					}	
					
				local cil = round(T1[`r3',`c'],.001)
				local ciu = round(T1[`r3',`c2'],.001)
				local cil = string(`cil', "%05.3f")
				local ciu = string(`ciu', "%05.3f")
				scalar conf_int = "[`cil' - `ciu']"
		
			local xcellcol: word `xc' of `excelcols'
			putexcel `xcellcol'`xr' = "`=pt_est'", hcenter
				local xr = `++xr'
				di `xr'
			putexcel `xcellcol'`xr' = "`=conf_int'", hcenter
				local xr = `++xr'
			}
		local xc = `xc'+1
		}
		
		putexcel B1:D1 = "Couples with No Previous Sons", merge border(bottom) hcenter
		putexcel F1:H1 = "Couples with At Least One Previous Son", merge border(bottom) hcenter
		putexcel B2 = "Pre-LLF", hcenter
		putexcel C2 = "Early LLF", hcenter
		putexcel D2 = "Late LLF", hcenter
		putexcel F2 = "Pre-LLF", hcenter
		putexcel G2 = "Early LLF", hcenter
		putexcel H2 = "Late LLF", hcenter
		putexcel A3 = "First Parity Births"
		putexcel A5 = "Second Parity Births"
		putexcel A7 = "Third + Parity Births"
		putexcel A2:D2 , border(bottom)
		putexcel F2:H2 , border(bottom)
		putexcel A8:H8 , border(bottom)
		
		putexcel save 
		
	********************************************
	* (3.2) Table 2
	********************************************
		putexcel set  `"$d_tab\Table_2.xlsx"'	, replace 
		local excelcols "B C D"
		
		local xc = 1
		forval c = 1(2)6 { 
			local xr = 4
			forval r = 1(3)6 {
				local c2 = `c'+1
				local r2 = `r'+1
				local r3 = `r'+2
				if T2[`r2',`c']<.1 {
					local lincom = round(T2[`r',`c'],.001)
					local lincom = string(`lincom', "%05.3f")
					scalar pt_est = "`lincom'**"
					}
				if T2[`r2',`c']<.05 {
					local lincom = round(T2[`r',`c'],.001)
					local lincom = string(`lincom', "%05.3f")
					scalar pt_est = "`lincom'**"
					}
				if T2[`r2',`c']<.001 {
					local lincom = round(T2[`r',`c'],.001)
					local lincom = string(`lincom', "%05.3f")
					scalar pt_est = "`lincom'***"
					}
				else if T2[`r2',`c']>=.1 {
					local lincom = round(T2[`r',`c'],.001)
					local lincom = string(`lincom', "%05.3f")
					scalar pt_est = "`lincom'"
					}	
					
				local cil = round(T2[`r3',`c'],.001)
				local ciu = round(T2[`r3',`c2'],.001)
				local cil = string(`cil', "%05.3f")
				local ciu = string(`ciu', "%05.3f")
				scalar conf_int = "[`cil' - `ciu']"
		
			local xcellcol: word `xc' of `excelcols'
			putexcel `xcellcol'`xr' = "`=pt_est'", hcenter
				local xr = `++xr'
				di `xr'
			putexcel `xcellcol'`xr' = "`=conf_int'", hcenter
				local xr = `++xr'
			}
		local xc = `xc'+1
		}
		
		putexcel B2:D3 = "-", hcenter
		putexcel B1 = "Pre-LLF", hcenter
		putexcel C1 = "Early LLF", hcenter
		putexcel D1 = "Late LLF", hcenter
		putexcel A2 = "First Parity Births"
		putexcel A4 = "Second Parity Births"
		putexcel A6 = "Third + Parity Births"
		putexcel A1:D1 , border(bottom)
		putexcel A7:D7 , border(bottom)
		
		putexcel save 
		
	********************************************
	* (3.3) Appendix Table A9
	********************************************
		mat TA9 = TA9_p1 \ TA9_p2
		putexcel set  `"$d_tab\Table_A9.xlsx"', replace 	
		local excelcols "B C D E F G"
		
		local xc = 1
		forval c = 1(2)12 { 
			local xr = 4
			forval r = 1(3)9 {
				local c2 = `c'+1
				local r2 = `r'+1
				local r3 = `r'+2
				if TA9[`r2',`c']<.1 {
					local lincom = round(TA9[`r',`c'], 0.001)
					local lincom = string(`lincom', "%05.3f")
					scalar pt_est = "`lincom'**"
					}
				if TA9[`r2',`c']<.05 {
					local lincom = round(TA9[`r',`c'], 0.001)
					local lincom = string(`lincom', "%05.3f")
					scalar pt_est = "`lincom'**"
					}
				if TA9[`r2',`c']<.001 {
					local lincom = round(TA9[`r',`c'], 0.001)
					local lincom = string(`lincom', "%05.3f")
					scalar pt_est = "`lincom'***"
					}
				else if TA9[`r2',`c']>=.1 {
					local lincom = round(TA9[`r',`c'], 0.001)
					local lincom = string(`lincom', "%05.3f")
					scalar pt_est = "`lincom'"
					}	
					
				local cil = round(TA9[`r3',`c'], 0.001)
				local ciu = round(TA9[`r3',`c2'], 0.001)
				local cil = string(`cil', "%05.3f")
				local ciu = string(`ciu', "%05.3f")
				scalar conf_int = "[`cil' - `ciu']"
		
			local xcellcol: word `xc' of `excelcols'
			putexcel `xcellcol'`xr' = "`=pt_est'", hcenter
				local xr = `++xr'
				di `xr'
			putexcel `xcellcol'`xr' = "`=conf_int'", hcenter
				local xr = `++xr'
			}
		local xc = `xc'+1
		}

		local xc = 1
		forval c = 1(2)12 { 
			local xr = 12
			forval r = 10(3)18 {
				local c2 = `c'+1
				local r2 = `r'+1
				local r3 = `r'+2
				if TA9[`r2',`c']<.1 {
					local lincom = round(TA9[`r',`c'], .001)
					local lincom = string(`lincom', "%05.3f")
					scalar pt_est = "`lincom'**"
					}
				if TA9[`r2',`c']<.05 {
					local lincom = round(TA9[`r',`c'], .001)
					local lincom = string(`lincom', "%05.3f")
					scalar pt_est = "`lincom'**"
					}
				if TA9[`r2',`c']<.001 {
					local lincom = round(TA9[`r',`c'], .001)
					local lincom = string(`lincom', "%05.3f")
					scalar pt_est = "`lincom'***"
					}
				else if TA9[`r2',`c']>=.1 {
					local lincom = round(TA9[`r',`c'], .001)
					local lincom = string(`lincom', "%05.3f")
					scalar pt_est = "`lincom'"
					}	
				local cil = round(TA9[`r3',`c'], .001)
				local ciu = round(TA9[`r3',`c2'], .001)
				local cil = string(`cil', "%05.3f")
				local ciu = string(`ciu', "%05.3f")
				scalar conf_int = "[`cil' - `ciu']"
		
			local xcellcol: word `xc' of `excelcols'
			putexcel `xcellcol'`xr' = "`=pt_est'", hcenter
				local xr = `++xr'
				di `xr'
			putexcel `xcellcol'`xr' = "`=conf_int'", hcenter
				local xr = `++xr'
			}
		local xc = `xc'+1
		}	
		
		putexcel B1:G1 = "Sex Ratio at", hcenter
		putexcel B2 = "Birth", hcenter
		putexcel C2 = "Age 1", hcenter
		putexcel D2 = "Age 2", hcenter
		putexcel E2 = "Age 3", hcenter
		putexcel F2 = "Age 4", hcenter
		putexcel G2 = "Age 5", hcenter
		putexcel A3 = "Second Births"
		putexcel A4 = "   Pre-LLF"
		putexcel A6 = "   Early LLF"
		putexcel A8 = "   Late LLF"
		putexcel A11 = "Third and Higher Order Births"
		putexcel A12 = "   Pre-LLF"
		putexcel A14 = "   Early LLF"
		putexcel A16 = "   Late LLF"
		putexcel A2:G2 , border(bottom)
		putexcel A17:G17 , border(bottom)
		
		putexcel save 
		
		********************************************
		* (3.4) Appendix Table A11
		********************************************
		lab var noson "No Previous Son"
		lab var period_2 "Early LLF Period"
		lab var period_3 "Late LLF Period"
		lab var period2_noson `"No Previous Son * Early LLF Period"'
		lab var period3_noson "No Previous Son * Late LLF Period"
		
		esttab preferred_1 robust_1 preferred_2 robust_2 preferred_3 robust_3 using `"$d_tab\Table_A11.html"' ///
		, replace cells(b(fmt(3)) ci(fmt(3) par("[" " - " "]"))) html label  alignment(center) ///
		keep(noson period_2 period_3 period2_noson period3_noson) eqlabels(none) nonumbers  ///
		order(noson period_2 period_3 period2_noson period3_noson)  ///
		mtitles("Preferred Specification" "Controlling for Family Size" "Preferred Specification" "Controlling for Family Size" "Preferred Specification" "Controlling for Family Size") collabels(none) ///
		mgroups("First Births" "Second Births" "Third Births", pattern(1 0 1 0 1 0) )
		
		
***********************************************
* (3.0) Compute population-level incidence 
*		of each sex selection strategy by 
*		sibship sex composition and LLF period
***********************************************
	preserve 
	
	* Set up some matrices to hold the results
	* Make sure no missing entries in stop_ or postnatal_ result sets
	foreach period in pre early late { 
		forval i = 1/3 { 
			forval j = 1/2 { 
			
				if postnatal_`period'[`i', `j']  < 0 { 
					mat postnatal_`period'[`i', `j']  = 0
					}				
				if postnatal_`period'[`i', `j']  == . { 
					mat postnatal_`period'[`i', `j']  = 0
					}
					
				if stop_`period'[`i', `j'] < 0 { 
					mat stop_`period'[`i', `j']  = 0
					} 
				if stop_`period'[`i', `j'] == .  { 
					mat stop_`period'[`i', `j']  = 0
					} 
					
				} 
			} 
		}
			
	************************************
	* (3.1) Parity and sex composition  
	* 		specific incidecne 
	************************************
	
	foreach period in pre early late { 
		mat incidence_any_cell_`period' = J(3, 2, . )
		mat incidence_stop_cell_`period' = J(3, 2, . )
		}
		
	foreach period in pre early late { 
		forval i = 1/3 { 
			forval j = 1/2 { 
				
				* Incidence of any sex selection
				mat incidence_any_cell_`period'[`i', `j'] = stop_`period'[`i', `j']*.515
				
				* Subtract postnatal selection for stopping rules use
				mat incidence_stop_cell_`period'[`i', `j'] 	= incidence_any_cell_`period'[`i', `j']-postnatal_`period'[`i', `j']
				
				}
			}	
		}
		
	* replace missing and negative with zeros
	foreach period in pre early late { 
		forval i = 1/3 { 
			forval j = 1/2 { 
				if incidence_any_cell_`period'[`i', `j'] == . { 
					mat incidence_any_cell_`period'[`i', `j'] = 0
					}
				if incidence_stop_cell_`period'[`i', `j'] == . { 
					mat incidence_stop_cell_`period'[`i', `j'] = 0
					}
				if incidence_stop_cell_`period'[`i', `j'] < 0 { 
					mat incidence_stop_cell_`period'[`i', `j'] = 0
					}
				}
			}
		}	

	************************************
	* (3.2) Parity specific incidence
	* 		- proportion of all births  
	************************************
	
	* Calculate the weights to apply to each group 
	* of births (% of parity-specific births occurring 
	* in each sibship sex composition group)
	
	foreach period in pre early late { 
		forval i = 1/3 { 
			scalar rowsum_`period'_`i' = 0 
			forval j = 1/2 { 
				scalar rowsum_`period'_`i' = rowsum_`period'_`i' + n_`period'[`i', `j']
				} 
			} 
	
		mat weight_`period' = J(3, 2, . )
	
		forval i = 1/3 { 
			forval j = 1/2 { 
				mat weight_`period'[`i', `j'] = n_`period'[`i', `j']/rowsum_`period'_`i'
				}
			}

		mat weighted_any_`period' = J(3, 2, . )
		mat weighted_stop_`period' = J(3, 2, . )
	
		forval i = 1/3 { 
			forval j = 1/2 { 
				mat weighted_any_`period'[`i', `j'] = incidence_any_cell_`period'[`i', `j']*weight_`period'[`i', `j']
				mat weighted_stop_`period'[`i', `j'] = incidence_stop_cell_`period'[`i', `j']*weight_`period'[`i', `j']
				}
			}
		}

	* Within-parity incidence 
	foreach period in pre early late { 
		forval i = 1/3 { 
	
			scalar incidence_any_row_`period'_`i' = 0 
			scalar incidence_stop_row_`period'_`i' = 0 

			forval j = 1/2 { 
		
				scalar incidence_any_row_`period'_`i' = incidence_any_row_`period'_`i' + weighted_any_`period'[`i', `j']
				scalar incidence_stop_row_`period'_`i' = incidence_stop_row_`period'_`i' + weighted_stop_`period'[`i', `j']
				} 
			} 		
				
		mat incidence_any_row_`period' = incidence_any_row_`period'_1 \ incidence_any_row_`period'_2 \ incidence_any_row_`period'_3 			
		mat incidence_stop_row_`period' = incidence_stop_row_`period'_1 \ incidence_stop_row_`period'_2 \ incidence_stop_row_`period'_3
		}	
	
	****************************************
	* (3.3) Population level incidence  
	* 		% of couples using each method
	****************************************
	
	* Use the parity specific incidence, 
	* weight for share of terminal births occurring 
	* at each parity in each pairty. 
	* Sum across all parities

	count if terminal == 1 & period == 1
		local count_pre = r(N)
	count if terminal == 1 & period == 2
		local count_early = r(N)
	count if terminal == 1 & period == 3
		local count_late = r(N)
	
	forval i = 1/2 { 
		count if terminal == 1 & order == `i' & period == 1
			local temp1 = r(N)
			local temp2 = `temp1'/`count_pre' 
			scalar order`i'_pre = `temp2'
		count if terminal == 1 & order == `i' & period == 2
			local temp1 = r(N)
			local temp2 = `temp1'/`count_early' 
			scalar order`i'_early = `temp2'
		count if terminal == 1 & order == `i' & period == 3
			local temp1 = r(N)
			local temp2 = `temp1'/`count_late' 
			scalar order`i'_late = `temp2'
		} 

		count if terminal == 1 & order >= 3 & period == 1
			local temp1 = r(N)
			local temp2 = `temp1'/`count_pre' 
			scalar order3_pre = `temp2'
		count if terminal == 1 & order >= 3 & period == 2
			local temp1 = r(N)
			local temp2 = `temp1'/`count_early' 
			scalar order3_early = `temp2'
		count if terminal == 1 & order >= 3 & period == 3
			local temp1 = r(N)
			local temp2 = `temp1'/`count_late' 
			scalar order3_late = `temp2'
		
		
		mat parity = 1 \ 2 \ 3 
		mat space = J(3,1,.)
		
		foreach period in pre early late { 
			mat pct_terminal_births_`period' = order1_`period' \ order2_`period' \ order3_`period' 
				
			mat incidence_any_`period' 	= parity, incidence_any_row_`period', pct_terminal_births_`period', space 
			mat incidence_stop_`period' = parity, incidence_stop_row_`period', pct_terminal_births_`period', space
		
			forval i = 1/3 { 
				mat incidence_any_`period'[`i', 4] 		=  incidence_any_`period'[`i', 2]*incidence_any_`period'[`i', 3]
				mat incidence_stop_`period'[`i', 4] 	=  incidence_stop_`period'[`i', 2]*incidence_stop_`period'[`i', 3]
				}
		
			scalar pop_incidence_any_`period' 	= 0 
			scalar pop_incidence_stop_`period' 	= 0 
		
			forval i = 1/3 { 
				scalar pop_incidence_any_`period' 	= pop_incidence_any_`period' + incidence_any_`period'[`i', 4] 	
				scalar pop_incidence_stop_`period' 	= pop_incidence_stop_`period' + incidence_stop_`period'[`i', 4] 	
				}	
			}
			
		mat population_incidence = pop_incidence_any_pre, pop_incidence_any_early, pop_incidence_any_late, pop_incidence_stop_pre, pop_incidence_stop_early, pop_incidence_stop_late 
		
	****************************************
	* (3.4) Population level incidecne  
	* 		% of births postnatally selected
	****************************************
	
	* Weight postnatal selection for share of births in each cell, 
	* weight parity-specific incidence for porportion of births in each parity 
	* Sum across all parities

	foreach period in pre early late { 
		forval i = 1/3 { 
			scalar rowsum_`period'_`i' = 0 
			forval j = 1/2 { 
				scalar rowsum_`period'_`i' = rowsum_`period'_`i' + n_`period'[`i', `j']
				} 
			} 

		mat weight_`period' = J(3, 2, . )
	
		forval i = 1/3 { 
			forval j = 1/2 { 
				mat weight_`period'[`i', `j'] = n_`period'[`i', `j']/rowsum_`period'_`i'
				}
			}				
	}
	
	foreach period in pre early late { 
		mat weighted_postnatal_`period' = J(3, 2, . )
		forval i = 1/3 { 
			forval j = 1/2 { 
				mat weighted_postnatal_`period'[`i', `j'] = postnatal_`period'[`i', `j']*weight_`period'[`i', `j']
				}
			}
		}
		
	foreach period in pre early late { 
		forval i = 1/3 { 
			scalar incidence_postnatal_row_`period'_`i' = weighted_postnatal_`period'[`i', 1] + weighted_postnatal_`period'[`i', 2] 
			} 			
		mat incidence_postnatal_row_`period' = incidence_postnatal_row_`period'_1 \ incidence_postnatal_row_`period'_2 \ incidence_postnatal_row_`period'_3 			
		}
		
	count if period == 1
		local count_pre = r(N)
	count if period == 2
		local count_early = r(N)
	count if period == 3
		local count_late = r(N)

	forval i = 1/2 { 
		count if order == `i' & period == 1 & terminal == 1
			local temp1 = r(N)
			local temp2 = `temp1'/`count_pre' 
			scalar order`i'_pre = `temp2'
		count if order == `i' & period == 2 & terminal == 1
			local temp1 = r(N)
			local temp2 = `temp1'/`count_early' 
			scalar order`i'_early = `temp2'
		count if order == `i' & period == 3 & terminal == 1
			local temp1 = r(N)
			local temp2 = `temp1'/`count_late' 
			scalar order`i'_late = `temp2'
		} 
	
		count if order >= 3 & period == 1 
			local temp1 = r(N)
			local temp2 = `temp1'/`count_pre' 
			scalar order3_pre = `temp2'
		count if order >= 3 & period == 2
			local temp1 = r(N)
			local temp2 = `temp1'/`count_early' 
			scalar order3_early = `temp2'
		count if order >= 3 & period == 3 
			local temp1 = r(N)
			local temp2 = `temp1'/`count_late' 
			scalar order3_late = `temp2'
			
	mat parity = 1 \ 2 \ 3 
	mat space = J(3,1,.)
	
	foreach period in pre early late { 
		mat pct_births_`period' = order1_`period' \ order2_`period' \ order3_`period' 
			
		mat incidence_postnatal_`period' = parity, incidence_postnatal_row_`period', pct_births_`period', space 
	
		forval i = 1/3 { 
			mat incidence_postnatal_`period'[`i', 4] 		=  incidence_postnatal_`period'[`i', 2]*incidence_postnatal_`period'[`i', 3]
			}
	
		scalar overall_incidence_post_`period' = 0 
	
		forval i = 1/3 { 
			scalar overall_incidence_post_`period' = overall_incidence_post_`period' + incidence_postnatal_`period'[`i', 4] 	
			}	
		}
		
	mat overall_incidence_postnatal = overall_incidence_post_pre, overall_incidence_post_early, overall_incidence_post_late 
	
	***************************************
	* (3.5) File clean up
	***************************************
	
		mat parity = 1 \ 2 \ 3 
	
		mat postnatal_pre = parity, postnatal_pre
		mat postnatal_early = parity, postnatal_early
		mat postnatal_late = parity, postnatal_late
	
		mat colnames postnatal_pre = "Parity" "NoSons" "Sons" 
		mat colnames postnatal_early = "Parity" "NoSons" "Sons" 
		mat colnames postnatal_late = "Parity" "NoSons" "Sons" 
		
		mat colnames incidence_postnatal_pre = "Parity" "WIParityIncidence" "Weight" "WeightedParIncidence"
		mat colnames incidence_postnatal_early = "Parity" "WIParityIncidence" "Weight" "WeightedParIncidence"
		mat colnames incidence_postnatal_late = "Parity" "WIParityIncidence" "Weight" "WeightedParIncidence"
	
	foreach period in pre early late { 
		mat incidence_any_cell_`period' = parity, incidence_any_cell_`period'
		mat incidence_stop_cell_`period' = parity, incidence_stop_cell_`period'

		mat colnames incidence_any_cell_`period' = "Parity" "NoSons" "Sons" 
		mat colnames incidence_stop_cell_`period' = "Parity" "NoSons" "Sons" 

		mat colnames incidence_any_`period' = "Parity" "WIParityIncidence" "Weight" "WeightedParIncidence" 
		mat colnames incidence_stop_`period' = "Parity" "WIParityIncidence" "Weight" "WeightedParIncidence"
		
		}
	
	mat colnames population_incidence = "AnyPre" "AnyEarly" "AnyLate" "StopPre" "StopEarly" "StopLate" 
	mat colnames overall_incidence_postnatal = "PostnatalPre" "PostnatalEarly" "PostnatalLate" 
	
	*************************************************************************************************************
	* (3.6) Output key results to data files
	*
	* - Any sex selective discontinuation incidence during each period 
	*		by sibship sex composition: Incidence_Any_bySons_pre (_early, _late)
	* - Sex selective stopping rules incidence during each period 
	*		by sibship sex composition: Incidence_Stop_bySons_pre (_early, _late)
	* - Postnatal selection incidence during each period 
	*		by sibship sex composition: Incidence_Postnatal_bySons_pre (_early, _late)	
	* - Overall incidence of any sex selective discontinuation (weighting for share 
	*		of births in each sibship sex composition group): OverallIncidence_Any_pre (_early, _late)
	* - Overall incidence sex selective stopping rules (weighting for share 
	*		of births in each sibship sex composition group): OverallIncidence_Stop_pre (_early, _late)
	* - Overall incidence of postnatal selection (weighting for share 
	*		of births in each sibship sex composition group): OverallIncidence_Postnatal_pre (_early, _late)
	*	
	* - Population-level incidence of any sex selective discontinuation and stopping rules per couple
	*		during each period: Population_Incidence_perCouple.dta
	* - Population-level incidence of postnatal selection per birth during each period: 
	*		Postnatal_Incidence_perBirth.dta
	*************************************************************************************************************
	
	cap mkdir `"$d_data\Sex Selection Incidence"'
	* These give the parity and period-specific incidence 
	foreach period in pre early late { 
		clear 
		svmat postnatal_`period', names(col) 
			save `"$d_data\Sex Selection Incidence\Incidence_Postnatal_bySons_`period'LLF.dta"', replace
		clear 
		svmat incidence_postnatal_`period', names(col) 
			save `"$d_data\Sex Selection Incidence\OverallIncidence_Postnatal_`period'LLF.dta"', replace
		}
	
	foreach period in pre early late { 
		clear
		svmat incidence_any_cell_`period', names(col)  
			save `"$d_data\Sex Selection Incidence\Incidence_Any_bySons_`period'LLF.dta"', replace
		clear
		svmat incidence_stop_cell_`period', names(col) 
			save `"$d_data\Sex Selection Incidence\Incidence_Stop_bySons_`period'LLF.dta"', replace
		clear
		svmat incidence_any_`period', names(col) 
			save `"$d_data\Sex Selection Incidence\OverallIncidence_Any_`period'LLF.dta"', replace
		clear
		svmat incidence_stop_`period', names(col) 
			save `"$d_data\Sex Selection Incidence\OverallIncidence_Stop_`period'LLF.dta"', replace
		} 	

		
	clear
	* This gives the population-level incidence of stopping rules and all discontinuation (percent of couples)
	svmat population_incidence, names(col)
		save `"$d_data\Sex Selection Incidence\Population_Incidence_perCouple.dta"', replace
				
	clear
	* This gives the share of births postnatally selected 
	svmat overall_incidence_postnatal, names(col)
		save `"$d_data\Sex Selection Incidence\Postnatal_Incidence_perBirth.dta"', replace 
	
	********************************************
	* (3.7) Table 3
	********************************************
	
	putexcel set  `"$d_tab\Table_3.xlsx"', replace 
	
		putexcel B1:D1 = "Couples with No Previous Sons", merge border(bottom) hcenter
		putexcel F1:H1 = "Couples with At Least One Previous Son", merge border(bottom) hcenter
		putexcel B2 = "Pre-LLF", hcenter
		putexcel C2 = "Early LLF", hcenter
		putexcel D2 = "Late LLF", hcenter
		putexcel F2 = "Pre-LLF", hcenter
		putexcel G2 = "Early LLF", hcenter
		putexcel H2 = "Late LLF", hcenter
		putexcel A3 = "First Parity Births"
		putexcel A4 = "Second Parity Births"
		putexcel A5 = "Third + Parity Births"
		
		forval i = 1/3 { 
			local xr = `i'+2
			local cell = round(incidence_stop_cell_pre[`i',2],.001)
				putexcel B`xr' = `cell', hcenter
			local cell = round(incidence_stop_cell_pre[`i',3],.001)
				putexcel F`xr' = `cell', hcenter
			local cell = round(incidence_stop_cell_early[`i',2],.001)
				putexcel C`xr' = `cell', hcenter
			local cell = round(incidence_stop_cell_early[`i',3],.001)
				putexcel G`xr' = `cell', hcenter		
			local cell = round(incidence_stop_cell_late[`i',2],.001)
				putexcel D`xr' = `cell', hcenter
			local cell = round(incidence_stop_cell_late[`i',3],.001)
				putexcel H`xr' = `cell', hcenter		
		}
		
		putexcel F3:H3 = "-", hcenter
		putexcel A2:D2 , border(bottom)
		putexcel F2:H2 , border(bottom)
		putexcel A5:H5 , border(bottom)
		
		putexcel save 
	
	********************************************
	* (3.8) Figure 4
	********************************************

	use `"$d_data\Sex Selection Incidence\Population_Incidence_perCouple.dta"', clear 
	
	rename AnyPre Any1
	rename AnyEarly Any2
	rename AnyLate Any3
	rename StopPre Stop1
	rename StopEarly Stop2
	rename StopLate Stop3
	
	gen n = 1 
	reshape long Any Stop, i(n) j(period)
	replace Any = Any*100
	replace Stop = Stop*100
	
	sum Stop if period==1 
		local stop_pre = r(mean)
		local stop_pre = round(`stop_pre', .01)
		local pre_pos = `stop_pre'-.3
	sum Stop if period==2 
		local stop_early = r(mean)
		local stop_early = round(`stop_early', .01)
		local early_pos = `stop_early'-.3
	sum Stop if period==3 
		local stop_late = r(mean)
		local stop_late = round(`stop_late', .01)
		local late_pos = `stop_late'-.3
	sum Any if period==3
		local any_late = r(mean)
		local postnatal = `any_late'-`stop_late'
		local postnatal = round(`postnatal', .01)
		local postnatal = string(`postnatal', "%04.2f")
		local postnatal_pos = `any_late'-.2 
		
	twoway (bar Any period, sort fcolor(black) barwidth(.75)) ///
	(bar Stop period, sort fcolor(gs10) barwidth(.75)) pcarrowi `postnatal_pos' 2.45 `postnatal_pos' 2.6, ///
	ytitle(Percent of Couples) xlabel(1 "Pre-LLF" 2 "Early LLF" 3 "Late LLF", noticks labsize(small)) ///
	ylabel(1 "1.0%" 2 "2.0%" 3 "3.0%" 4 "4.0%" 5 "5.0%" 6 "6.0%" 7 "7.0%" 8 "8.0%", grid glcolor(black) labsize(small)) ///
	xtitle("") legend(pos(6) rows(1) order(2 "Stopping Rule Use" 1 "Postnatal Neglect" ) size(small)) ///
	text(`pre_pos' 1 "{bf: `stop_pre'%}" `early_pos' 2 "{bf: `stop_early'%}" `late_pos' ///
	3 "{bf: `stop_late'%}" `postnatal_pos' 2.35 "{bf: `postnatal'%}" , size(small))
	
	graph export `"$d_fig\Fig_4_SexSelection.jpg"', replace
	